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^ ■ Heuristic methods for solution of problems in the NP-Complete class of 

decision problems often reach exact solutions, but fail badly at "phase bound- 



j^ \ aries," across which the decision to be reached changes from almost always 

having one value to almost always having a different value. We report an 
^ . analytic solution and experimental investigations of the phase transition that 

o, 

O , occurs in the limit of very large problems in K-S AT. The nature of its " random 

first-order" phase transition, seen at values of K large enough to make the 

(•^ ■ computational cost of solving typical instances increase exponentially with 

OO . 

f^ , problem size, suggests a mechanism for the cost increase. There has been 

o 

evidence for features like the "backbone" of frozen inputs which characterizes 
^* ■ the UNSAT phase in K-SAT in the study of models of disordered materials. 



a 



o 



but this feature and this transition are uniquely accessible to analysis in K- 
SAT. The random first-order transition combines properties of the 1st order 
(discontinuous onset of order) and 2nd order (with power law scaling, e.g. of 
the width of the critical region in a finite system) transitions known in the 



k> i physics of pure solids. Such transitions should occur in other combinatoric 

C^ \ problems in the large N limit. Finally, improved search heuristics may be 

developed when a "backbone" is known to exist. 



Constraint satisfaction, the automated search for a configuration of a complex system 
which satisfies a set of rules or inequalities, is often difficult, and occurs widely in practice. 
The simplest example of a constraint satisfaction problem, K-SAT |I|,^ , is commonly used 
as a testbed for heuristic algorithms intended for wider use and was the first problem proved 
to be in the complexity class NP-Complete |]^J§], in which the worst case instances are 
believed to always require computing effort exponential in A^, the number of input degrees 
of freedom. In random K-SAT, the system parameters are a string of A^ bits, and the rules 
to be satisfied are a set of M clauses. If the string consists of bits {xi = 0, l}i=i,...,Af, we 
construct an instance of K-SAT by first randomly choosing K distinct possible indices i and 
then, for each of them, a literal Zi (i.e. the corresponding Xi or its negation Xj with equal 
probability). A clause C is defined as the logical OR of the K literals chosen. Next, we 
repeat this process to obtain M independently chosen clauses {Ce}e=i,...,M and ask for all of 
them to be true at the same time, i.e. we take the logical AND of the M clauses. This gives 
a formula F, which may be written as 

M M / K \ 

F=Ac, = /\ lyzn , (1) 

where A and V stand for the logical AND and OR operations respectively. An assignment of 
the {xj}'s satisfying all clauses is a solution of the K-SAT problem. If no such assignment 
exists, F is unsatisfiable. The formulae F constructed at random keeping the ratio a = M/N 
constant as M, N ^ oo provide a natural ensemble of test problems, with a characterizing 
whether the F are typically under- or over-constrained. 

The value of K is important. 2-SAT can be solved by a linear time algorithm [0. There 
is a critical value of a, ac{2) = 1, below which the likelihood of an F being UNSAT vanishes 
in the limit A^ -^ oo, and above which it goes to 1. For A' > 3, K-SAT is NP-Complete 
and rigorous results are few. Computer experiments P,^ on K-SAT for K = 2,3 and higher 
have located the phase transition and provided critical exponents for the sharpening of the 
critical region which occurs with increasing sample size A^. 

The technique used, finite-size scaling, will be discussed in more detail below. It has 
recently been proven [|l^ that certain monotonic properties (such as SAT) in combinatorial 
ensembles do have sharp thresholds. "Sharp" means in our case that for any a < a^K, N), 
the probability that a formula in K-SAT can be satisfied goes to 1 as A^ — * oo , while for any 



a > ac{K,N) this probability tends to 0. While Friedgut's result leaves open the question 
of whether a^K, N) has a limit as A^ ^ oo, experiment suggests that this is the case for 
K-SAT. 

One reason for the recent interest in the threshold is the growing recognition that "easy- 
hard-easy" patterns of computational cost are characteristic of heuristic algorithms devel- 
oped to solve problems whose cost in the worst case increases exponentially with problem 
size N, and that the hardest instances occur near |Tl| phase boundaries like Oc |T^-|14[. 

There is a strong analogy between these problems and the properties of disordered mate- 
rials, alloys or even glasses, studied by constructing models whose energy function expresses 
the constraints |T5[. Strongly disordered models with conflicting interactions similar to the 



randomly negated literals in K-SAT are known as "spin glasses" |T^. Fu and Anderson |0] 
first conjectured that spin glasses are the models underlying NP-Complete decision problems 
and vice-versa. They cite the example of weighted graph partitioning, which is equivalent 
to an Ising or Potts spin glass. 

A technique of calculating expectation values of observables in random many-parameter 



systems, called the replica method fT^ predicts that ordering of a new type is possible in the 
presence of microscopic randomness. In the replica method, the calculation of the average 
over the disorder leads to an effective energy or cost function which describes many identical 
copies of one instance of the model system with the dynamical variables (usually Ising spins 
taking on values +1 or -1) in different replicas coupled by some non-linear function. The 
onset of ordering can be identified by the fact that a single stable state occurs in multiple 
replicas. This corresponds to a physical ground state that is highly irregular in structure 
due to the randomness of the problem. A more subtle kind of ordering (called Replica 
Symmetry-Breaking) occurs when distinct stable states of this sort are found in different 
subsets of the replicas, signalling that there may be infinitely many distinct ground states 
with energies infinitesimally close to the true ground state. The extent of or absence of 
this new kind of order can be quantified by an "order parameter," which in general emerges 
from the replica formulation. We describe both types of result in sections I and II. Details 
of the calculations for K-SAT are given in section II. While the replica procedures are not 
rigorous, certain steps can be proven, and some results have been verified by other means. 



We discuss these issues in section II. 

Spin glass models with realistic connectivity are difficult to explore experimentally, either 
on real substances or by computer simulation of models in thermal equilibrium at some finite 
temperature. Experimental study of the "easy-hard-easy" phase transitions in combinatorics 
is more tractable. Although these are spin glass models, the properties of interest are 
ground state properties, and a large body of model-specific heuristics exists, which gives 
powerful means of exploring these ground states. We have previously applied replica methods 
and determined characteristics of the 3-SAT transition [|I^-|2D|. In section III, we report 
additional results which provide new insights, hopefully of use to both fields. Note, however, 
that the methods of statistical physics predict the most probable, or typical, behavior of a 
system with many degrees of freedom, so we shall be describing the typical complexity of 
K-SAT, not its worst case. 

I. MIXTURES OF if = 2 AND 3: OVERVIEW OF RESULTS 



In order to understand what occurs between K = 2 and K = 3, we have studied \TL 
formulae containing mixtures of 2- and 3-clauses: consider a random formula with M clauses, 
of which (1 — p)M contain two literals and pM contain 3 literals, with < p < 1. This 
"2 +p-SAT" model smoothly interpolates between 2-SAT {p = 0) and 3-SAT {p = 1). The 
problem is NP-complete, since any instance of the model for p > contains a sub-formula 
of 3-clauses. But our interest here is in the complexity of "typical" problem instances. 

We seek ac(2 -|- p), the threshold ratio M/N of the above model at fixed p. We know 
ac(2) = 1 and adS) ~ 4.27. F cannot be almost always satisfied if the number of 2-clauses 
(respectively 3-clauses) exceeds A^ (resp. q;c(3)A^). As a consequence, the critical ratio must 
be bounded by ac{2 + p) < min ij^, ^^^j. 

The 2 -I-J9-SAT model can be mapped onto a diluted spin glass model with A^ spins Sf. 
Si = 1 ii the Boolean variable Xi is true. Si = —1 ii Xi is false. Then, to any configuration is 
associated an energy E, or cost-function, equal to the number of clauses violated. Random 
couplings between the spins are induced by the clauses. The most important result of the 
replica approach ^T^J is the emergence, in the large N, M limit and at fixed p and a, of order 
parameters describing the statistics of optimal assignments, which minimize the number of 



violated clauses. In this section, we give an overview of the resuhs from statistical mechanics. 
The next section gives a more detailed description of the analysis. 

Consider an instance of the 2 + p-SAT problem. We use the N'gs ground state configu- 
rations to define 



^. = ^^ E s! (2) 

9= 
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the average value of spin Si over all optimal configurations. Clearly, rrii ranges from —1 to +1 
and rrii = —1 (respectively +1) means that the corresponding Boolean variable Xj is always 
false (resp. true) in all ground states. The distribution P{m) of all rrii gives the microscopic 
structure of the ground states. The accumulation of magnetizations m around ±1 represents 
a "backbone" of almost completely constrained variables, whose logical values cannot vary 
from solution to solution, while the center of the distribution P{m ~ 0) describes weakly 
constrained variables. The threshold ac will coincide with the appearance of an extensive 
backbone density of fully constrained variables Xj, with a finite probability weight at m = ±1. 
A simple argument shows that the backbone must vanish when a < ac- Consider adding 
one clause to a SAT formula found below ac- If there is a finite fraction of backbone spins, 
there will be a finite probability that the added clause creates an UNSAT formula, which 
cannot occur. 

For a < ac, the solution exhibits a simple symmetry property, usually referred to as 
Replica Symmetry (RS), which leads to an order parameter which is precisely the magneti- 
zation distribution P{m) defined above. An essential qualitative difference between 2-SAT 
and 3-SAT is the way the order parameter P{m) changes at the threshold. This discrepancy 
can be seen in the fraction f{K, a) of Boolean variables which become fully constrained, at 
and above the threshold. As said above, f{K, a) is identically null below the threshold. For 
2-SAT, f{2,a) becomes strictly positive above etc = 1 and is continuous at the transition : 
/(2, 1^) = /(2, 1+) = 0. On the contrary, /(3,a) displays a discontinuous change at the 
threshold : fi3,a~) = and /(3,a+) = /e(3) > 0. 

While for the continuous transitions, the exact value of the threshold can be derived 
within the RS scheme, for the discontinuous case the RS prediction gives only upper bounds. 
The exact value of the threshold can be predicted only by a proper choice of the order 
parameter at the transition point, i.e. by a more general symmetry breaking scheme, a 



problem which is still open. However, the predictions of the RS equations, such as the 
number of solutions, remain valid up to a^-, and the RS prediction for the nature of the 
threshold should be qualitatively correct. 

For the mixed 2 + p-SAT model, the key issue is therefore to understand how a discon- 
tinuous 3-SAT-like transition may appear when increasing p from zero up to one and how 
it affects the computational cost of finding solutions near threshold. Applying the method 
of ref. |T9|, we find for p < po {po = 0.41), there is a continuous SAT/UNSAT transition at 
ac{2+p) = j^ . This has recently been verified by rigorous analysis up to p = 0.4 |^2[. The 
RS theory appears to be correct for a < ac{2+p), and thus gives both the critical ratio and 
the typical number of solutions, as in the K = 2 case. The SAT/UNSAT transition should 
coincide with a replica symmetry breaking transition, as discussed in [0. So, for p < Po, 
the model shares the characteristics of random 2-SAT. 

For p > Po, the transition becomes discontinuous and the RS transition gives an upper 
bound for the true ac{2 +p). The RS theory correctly predicts a discontinuous appearance 
of a finite fraction of fully constrained variables which jumps from to fc when crossing the 
threshold ac{2 +p). However, both values of /c(2 + p) and Oc are slightly overestimated, 
e.g. for p = 1, a^^{3) ~ 4.60 and /c^'^(3) ~ 0.6 whereas experiments give ac(3) — 4.27 
and /c(3) ~ 0.5. A replica symmetry breaking theory will be necessary to predict these 
quantities. For p > pq, the random 2+p-SAT problem shares the characteristics of random 
3-SAT. 

This transition differs from phase transitions in most ordered materials in that the ground 
state is highly degenerate at Oc- The entropy, i.e. the logarithm (base 2) of the typical 
number of optimal solutions, can be computed exactly within the RS scheme for any p and 
a < oiJ(l + p). The entropy at the transition point decreases as a function of p, from 0.56 
for p = to 0.13 for p = 1, as plotted in Fig. 1. 

II. STATISTICAL MECHANICS ANALYSIS OF THE 2+P SAT MODEL 

In this section we describe the analytical calculation of the typical ground state properties 
of the 2+P-SAT model using the replica method. A brief discussion concerning rigorous 
results and prospects for making the replica results rigorous is also included. 



A. The energy-cost function 

The logical variables Xj can be represented by A^ binary variables Si, called spins, through 
the one-to-one mapping Si = —1 (respectively +1) if Xj is false (resp. true). We then encode 
the random clauses into a M x N matrix Cn in the following way : Cu = — 1 (respectively 
+ 1) if the clause Q includes Xi (resp. Xj), Cu = otherwise. Note that Y^iLi CuSi gives the 
the net number of literals satisfying clause £. Consider now the cost-function i5[C, S] defined 
as the number of clauses that are not satisfied by the logical assignment corresponding to 
configuration S. 
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where S[.; .] denotes the Kronecker function, which is 1 if its arguments are equal, zero 
otherwise. The minimum (or ground state) E[C] of E[C, S], the lowest number of violated 
clauses that can be achieved by the best possible logical assignment [|1^ , is a random variable 
which becomes totally concentrated around its mean value <^ E[C] ^ in the large size 
limit [^ . The latter is accessible through the knowledge of the averaged logarithm of the 
generating function 

Z[C] = Y.ewi-E[C,S]/T) (4) 

s 

since <^ E[C] ^= —T <^ logZ[C] ^ +0{T'^) when the auxiliary parameter T is eventually 
sent to zero. Since -C E[C] ^= in the SAT region and is positive in the UNSAT phase, 
calculating <^ E[C] ^ locates ac{K). 

B. The average over the disorder 

The calculation of the average value of the logarithm of Z from Eq. (^ is an awkward 
one. To circumvent this difficulty, we compute the n*^ moment of Z for integer- valued n and 
perform an analytical continuation to real n to exploit the identity ^ ^[0]" ^= 1 + n ^ 
logZ[C] ^ +0{n'^). The n^^ moment of Z is obtained by replicating n times the sum over 
the spin configuration S and averaging over the clause distribution |19| 



<Z[C]">= Y. <exp(-5]E[C,S"]/T) > 

Si.S2....,S" 



(5) 
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The average over the clauses can be performed because their probabihty distributions are 
uncorrelated. We obtain 



ai\il-p)M 



« z[cr »= E (C2n) 

Si,S2,...,S" 

where each factor is defined by {K = 2,3) 
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We stress that <^ . ^ now denotes the unbiased average over the set of S'^f^j vectors of A^ 
components C, = 0, ±1 and of squared norm equal to K. 
Resorting to the identity, 



N 



n s[st^-Ci_ 



6 J2C^St;-K 
.1=1 

one may carry out the average over disorder in eq.(^ to obtain 

1 1 ^ f 1 " -^ 1 

^ Ci,...,Ck=±1 Ji,...,iK=l l> a=l£=l ) 



(8) 



(9) 



to the largest order in A^. 

It is crucial to remark that C-ft:[5''^] in (^) depends on the nxN spins only through some 2" 
quantities c{a) labelled by the vectors a with n binary components; c{a) equals the number 
(divided by A^) of labels i such that S^ = a", Va = 1, . . . ,n |]2^. Indeed, one can rewrite 
a[51 = Ck[c\ with 

Ck[c] = ^ E E ci-C^a^)...ci-CKffK)exp\-^J2l[6K,l]] . (10) 

Notice that c(cr) = c(— a) due to the even distribution of the disorder C. 
Introducing the effective energy function. 
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we may rewrite the n^^ moment of the generating function Z (^ as 

«z"»=/nrfc(a)e-^^//w E nM^(^)-^^n^[5r;^i • (12) 

S Si,S2,...,S" a \ ^^ i=la=l J 

The sum over the spins in the last term of the above equation can be computed, and gives 
rise to a combinatorial factor 

m 



expf-iVEc(a)lnc((T)'| , (13) 



n.-(A^c(c?))! 

to the leading exponential order in A^. As a consequence, the ^^ moment of Z using the 
Laplace method is <^ Z" ^~ exp(A^ Fmax) where Fmax is the maximum over all possible cs 



of the functional 



F[c] = - E c(a) logc(a) - TE,^^\c\ , (14) 



with the constraint 



Ec(a) = l. (15) 



C. The replica symmetric theory 

The optimisation conditions over F\c\ provide 1^ coupled equations for the cs. Notice 
that F is a symmetric functional, invariant under any permutation of the replicas a. A 
maximum may thus be sought in the so-called replica symmetric (RS) subspace of dimension 
n + 1 where c(cr) is left unchanged under the action of the symmetric group. Within the 
RS subspace, the occupation fractions may be conveniently expressed as the moments of a 
probability distribution P{m) over the range — 1 < m < 1 ||19|| . 

/■I " /l + mcr"\ 
c(ai,a2,...,cr„) = / dmP^m)]^!^ ) . (16) 

P{m) is the distribution of Boolean magnetizations previously introduced in the paper. 

At this stage of the analysis it is possible to perform the analytic continuation n — *> 0, 
since all the functionals have been expressed in term of the generic number of replicas n. 
Such a process lead to a self-consistent functional equation for the order parameter P{m). 
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In the limit of interest T ^ 0, in order to properly describe the accumulation of the Boolean 
magnetization to the border of its domain {m G [—1,1]), it is convenient to introduce the 
rescaled variables z, implicitly defined by the relation m = tanh(z/T). Calling R{z) the 
probability distribution of the zs, we obtain 

/OO (j^y^ f POO 

-— cos('uz) exp < — «(1 — p) + 2«(1 — p) / c/2;i_R(2;i) cos('umin(l, 2:1)) 
-00 zvr I JO 

3 /""^ 1 

— -ap + 3ap dzidz2R{zi)R{z2) cos{uram(l, zi, Z2))> . (17) 
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As discussed in detail in ref. |T9|, the above type of equations admit an infinite sequence 



of rapidly converging exact solutions of the form 

R{z) = f: n6(z-^] . (18) 

£=-00 \ y/ 

In the above expression, 1/q is the resolution of the rescaled variable z which eventually goes 

to zero. Equation (|1^, leads to the following set of coupled equations for the coefficients 

r/s 



re 



J^ - cosm exp Y: 7,(cos(j0) - 1) (19) 



for all £ = 0, . . . , g — 1 where 

7j/a = 2 (1 - p) Tj + 3 p Tj 1 - ro - 2 ^ r^ - rJ , Vj = 1, . . . , g - 1 

^Ja = {l-p)ll-ro-2j2rA+^p fl-ro-2'^rJ . (20) 



By looking for the value of a at which the internal energy (|TT]) becomes positive, we are 
able to recover the results discussed previously in the text. For p < po, the transition is 
continuous and the solution of the equations up to ac is simply tq = 1, r^ = (/ = 1, ..., g). 
At ac, the coefficients r^ become continuously positive. For p > pq, the coefficients r^ jump 
discontinuously to a finite value beyond ac- It follows that in order to find the point where 
the discontinuous transition first takes place, one should look, within the RS scheme, for the 
point po at which the derivative of the order parameters re at ac diverge. 

We may expand the saddle point equations (0,^) to the second order in the parameters 
ri and s = 1 — tq. We find 
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(22) 



a^(l-p)^ 



q-l 



E^I + (2-E^^) 



^«^(i-p)V . 



Li=i «=i 

The analysis of the hnear terms in the above equations shows that the threshold is given by 

1 



ac{2+p) 
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iO<p<Po). 



(23) 



Next, we expand around the latter by posing a = -^ — \- x, r^ = B^x and s = Ax. At 
the critical point po, the above quantities {Be, A} should diverge in order to have a first 
order jump when x -^ 0^. We then assume that Be = XeA, with Xe = 0(1) and A -^ cxd, 
discarding irrelevant O(x^) corrections to the order parameters. We find q equations for Pq 
and A^, £ = 1, . . . , g — 1. 

2 



and, for 
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Though we have not been able to find an exact solution to (0,^) , the above equations can 
be easily solved iteratively, leading to a value of po — 0.41. A more detailed discussion of 



the equations for po is given in [^ , in which the connection with the calculation of ref . ^2 
is made explicit by showing that po = 2/5 is indeed a lower bound for pq. 

The exactness of the above results depends on the validity of the RS assumption in- 
troduced in the functional saddle-point equations. For p < pq, such an assumption turns 



10 



out to be correct, leading to a threshold value which as been proven to be exact also by 
other methods P2| For p > Pq, the change in the order parameter P{m) (or R{z)) at the 



threshold becomes discontinuous and the solutions of the RS equations account only for an 
upper bound of the true threshold. The exact solution of the self-consistency equations lies 
outside the RS subspace to which we have restricted our analysis. The exact determination 
of the SAT/UNSAT threshold for discontinuous transitions requires the introduction of a 
more general (and much more complicated) symmetry breaking scheme in the equations, 
the so called Replica Symmetry Breaking, which embodies the RS subspace as a particular 
case. 

It is worth noting that in the SAT region a < ac(p); the RS theory is believed to be 
exact and allows for the estimation of quantities of interest such as the typical number of 
solutions or the probability distribution of the variables over all solutions. Some rigorous 
probabilistic results are given in 



D. Comments on the replica approach 

In the previous paragraphs, replicas are introduced as a trick to compute the average 
value ^ logZ[C] ^ from the integer moments ^ Z[C]" ^ of the generating function (4). 
As long as the number of variables N is finite, a{n) =<^ Z[C]" ^ is an analytic function of n 
and grows less than exponentially at large n, a{n) < (2^)". Invoking a theorem of Carlson, 
a{n) is uniquely known from its values at the nonnegative integers. Therefore the analytic 
continuation to real n — > is unambiguous. However, due to the saddle-point calculation 
of Section II. B the limit A^ ^ oo is made first and the analytic continuation requires the 
introduction of some additional hypotheses. 

The most natural continuation scheme, called Replica Symmetry (RS) has recently been 



shown to be exact at high temperature T [P6| , p7| for the K-SAT model. Though not explicitly 



proven in p6| , p7| , it is reasonable to think that RS should also hold at T = for small ratios 



of clause per variable a. Indeed, in a simple constraint satisfaction model, the RS hypothesis 



has been shown |28] to be exact in the range < a < ac giving thus access to the exact 



value of the threshold Oc p9[ . Note that self-consistency criteria of the calculation of the 



local stability of the saddle-point over c(a) found in Section II. B are satisfied by the RS 
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hypothesis in this range. 

What happens when RS fails, e.g. above (respectively in the vicinity of) ac for 2- 
SAT (resp. 3-SAT) ? Among all models for which RS fails, the so called Random Energy 
Model (REM) p^] has been the only one rigorously solved so far. Its exact solution can 
be reproduced |^ using another Replica Symmetry Broken (RSB) scheme designed by 
Parisi fl^. Within such a scheme, the analytic continuation is performed by an iterative 
hierarchical procedure characterized by a closed algebraic structure at each stage of the 
hierarchy |]T^. The known random mean-field models (i.e. models with a complete graph 
of interaction) appear to be divided in two main classes. A first one for which the complete 
solution requires an infinite iteration of the Parisi scheme (e.g. the SK model [0]) and 
a second one for which the first step already provides the correct result. In the latter 
case, the successive steps of the RSB scheme lead to saddle point equations having as 



solutions the first step result |^. It is worth noting here that in the random 3-SAT problem, 
like several random mean-field models known to exhibit a discontinuous transition and be 
solvable by the one step RSB ansatz, the energy can be expressed as a sum of products of 
up to three spins S^. Therefore, the use of the one-step RSB hypothesis has promise for 
analyzing the SAT/UNSAT transition of 3-SAT. We expect that there will be differences 
resulting from the fact that the 3-SAT model involves a sparse graph of interactions. This 
permits heterogeneous orderings not possible in the mean-field models, for which all degrees 
of freedom are frozen in the ordered phase. 

III. EXPERIMENTAL TESTS 

We have run experiments to test the theory between K = 2 and K = 3, finding thresholds 
and assigning an exponent u for the narrowing of the critical region by finite-size scaling. 
The data obtained for this study is collected and presented in Figs 2a and 2b, which show the 
fraction of formulae that are unsatisfiable as a function of a for sample sizes from A^ = 50 to 
the largest practical size for each value of p. To obtain the data in Figs 2a and 2b, we take 
a sample of formulae with the desired value of p, and for each formula, starting well inside 
the SAT phase, add clauses until the formula first becomes unsatisfiable. The cumulative 
distribution of the values of a at which this occurs provides the curves in Figs 2. From 
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10,000 to 40,000 formulae were studied for each value of A^ and p shown. For cases with 
p > Po, our scripts used as their core the TABLEAU implementation of the Davis-Putnam 
search algorithm [^,Q. For cases with p < po a variant called MODOC was used ||3^. This 



adds binary resolution to eliminate 2-clauses by the relation (p V g) A {{~'q) V r) = pM r. 

Finite-size scaling of the critical region is done by plotting all quantities against the 
rescaled variable 

y = NHa- a,{K, N))/a,{K, N), (26) 

which "opens up" the critical region in samples with large N so that data from all sizes 
collapse onto a single universal curve [Q. If a^K, N) were constant, all the curves for a given 
value of K would pivot about a single point, adK, oo). This occurs in the finite-size scaling 
analysis of the random graph ensemble and is a good approximation at large K for K-SAT 
0. But there are significant additional size dependences present for small values of K, as 
evidenced in Fig. 3, which shows K = 2, {p = 0.0) on a fine scale. The successive crossings 
of pairs of curves for increasing values of N provides a rough measure of ac{K, N) (e.g. 
estimate a^K, 50) as the point where the fraction UNSAT for A^ = 50 crosses the fraction 
UNSAT for A^ = 100). If we make the ad hoc assumption that the added size dependence 
is due to the variation with A^ of ac{K,N), then we do not have to specify ac{K,N) for 
each A^. The data reduction required is to choose values of a^K, oo),!/ for which all the 
transformed curves are parallel, that is, shifted by {ac{K,N) — a^K, oo))/ a c{K,oo). An 
example of such a reduction is shown in Fig 4, for the case p = 0.0. Using this methodology, 
we obtained rescaled curves for all of the data which varied smoothly with p, as shown in 
Fig. 5. 

We find (Fig. 6) good agreement between the observed and predicted values of Oc, with 
an error which increases slowly from p = 0.41 to p = 1.0. We also show the bounds obtained 
by rigorous methods in Fig. 6. Lower bounds are obtained by showing that some analyzable 
algorithms, such as unit clause propagation [0 find SAT solutions with a finite probability 
|53|. Upper bounds use the fact that the probability of finding a satisfying assignment is 



bounded by the expected number of solutions. Refined versions of this argument |36,37 



partially eliminate the high degeneracy of some solutions [^]. Both methods have been 



applied to (2 +p)-SAT in p2[, yielding the dashed and dotted lines plotted in Fig. 6. 
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In figure 7, we show values of v obtained from finite-size scaling analysis described above. 
Below po; the exponent v is roughly constant and equal to 2.8, the value found for 2SAT. 
This indicates that the critical behavior along the second order transition line in Fig. 6 
is dominated by the 2-clauses in the formulae. When we include additional corrections to 
scaling in y{a,N) and in the probability of UNSAT, following the classic prescription [Q, 
we find that v may be as large as 3, the value that occurs in the percolation transition for 
undirected random graphs ||39|. The UNSAT phase for K = 2 is one in which "constraint 



loops" become so ubiquitous that almost certainly there is some literal that implies its 
converse. It is likely that the 2SAT transition results from percolation of these loops, and is 
in the same universality class as random graph percolation, differing only in the corrections 
to its scaling behavior. 

Above po, T-' drops rapidly to 1.5. For 7^ ^ 3, the values of z/ tend to 1.0, a result which 



can be understood in the "annealed" limit discussed in [IS, 3^. 

It is surprising that finite size scaling holds in the presence of discontinuous behavior of 
the order parameter characterizing the backbone. But this discontinuity is accompanied by 
smooth behavior of other thermodynamic quantities, e.g. entropy, as first discussed by Gross 
and Mezard in [^ . First order transitions in pure solids involve two (or a finite number of) 



phases and do not exhibit critical fluctuations or scaling laws with non-integer exponents. 
The random first-order transition taking place for p > po, into an infinite number of distinct 
ground states, displays features of both first and second orders. This mixed behaviour has 



also recently been observed in random-field models |^ . 

Previous work showed that the cost of running the best heuristics, depth-first search 
algorithms with backtracking [^ increased exponentially with A^ for any value of a for 



K = 3, with a prefactor that could be mapped into a universal function by plotting it as 
a function of y [|1^. The cost was maximized at a,5{K,N), so we have obtained cost data 
at this value of a for p = 0., .2, .4 and .6 over the range of N that could be searched. The 
plot in Fig. 8a shows that the median cost increases linearly with N for p < pq. It increases 
dramatically over a smaller range of A^ for p > pq. Fig. 8b confirms that this increase is 
exponential already for p = 0.6. 

Discontinuous nucleation of UNSAT regions due to the breakdown of replica symmetry 
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and the "backbone" of frozen spins conveniently explain the apparent inevitably high cost 
of heuristic search near the phase boundary. Heuristics which proceed by "asserting" the 
possible value of a spin make early mistakes by mis-asserting a backbone spin, and take 
a long time to backtrack to correct their mistakes. Even if the backbone can be identified 
before the depth-first search begins, the problem that remains is one of organizing the search 
over the remaining spins which lie on the boundaries of "nuclei" or partial solutions to find 
the lowest energy arrangements of the whole solution, also involving much wasted effort to 
explore an exponential subspace. 

The experiments, shown in Figs. 9a {K = 2) and 9b {K = 3), confirm that the appear- 
ance of the backbone is discontinuous for K = 3, and support the prediction of a continuous 
appearance of the backbone for K = 2. Above the threshold, the fraction of frozen spins 
found in small samples by exhaustive enumeration to locate all ground states is relatively 
insensitive to N. At and below the threshold, the fraction of frozen spins decreases rapidly 
with increasing N. While the samples which could be studied are too small to permit ex- 
trapolation, the results are consistent with < f{K, a) > vanishing below ac- 

The existence of a "backbone" has previously been reported in the traveling salesman 
problem [^ , with only a few bonds differing in many near-optimal tours. This observation 



has recently been turned to advantage by heuristics |^2[ which identify the backbone links 
and concentrate attention on the small subproblems which remain. This may prove to 
be a generally valid approach. Efficient means of finding the backbone will be specific to 
each problem type, but should nonetheless provide a step ahead in algorithm efficiency. 
Moreover, many worst-case NP-complete problems occurring in practice contain a mix of 
tractable and intractable constraints. Our results suggest that search algorithms that exploit 
as much of the tractable structure as possible may in fact scale polynomially in the typical 
case. In much of the work on search methods in, for example. Artificial Intelligence and 
Operations Research, one already informally follows the methodology of exploiting tractable 
problem structure in worst-case intractable problems. However, our hybrid model provides 
the first formal explanation why such a methodology can work so well in practice. Below 
a certain threshold fraction of intractable constraints, the overall behavior is dominated by 
the tractable substructure of the problem, leading to an overall efficient, polynomial time 
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solution method. 
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1. Ground state entropy at adp) versus p, predicted by the RS theory of PU| , PU[ . For 
p < po, adp) = 1/(1 — p). For p > pq, we have used the estimates of adp) obtained by 
finite-size scahng. 



p = 0.0, 0.2, 0.4 
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2a. Raw data used in this study for p = 0.0, 0.2, and 0.4. Vertical hues mark the 
thresholds for each value of p. For p = 0.0{K = 2), data are plotted for N = 50, 100, 200, 
500, 1000, 2000, 5000, 7500, and 10000. For p = 0.2, values of A^ are 100, 200, 500, 1000, 
2000, 5000, and 7500. For p = 0.4, values of A^ are 100, 200, 500, 1000, 2000, 3500, and 
5000. 
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p = 0.5, 0.6,0.8, 1.0 
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2b. Raw data for p = 0.5, 0.6, 0.8 and 1.0. Thresholds marked are determined from the 
RS theory (an overestimate). For p = 0.5, values of A^ are 50, 100, 150, 250, 500, 1000, 
1500, 2000, and 2500. For p = 0.6, values of A^ are 50, 100, 150, 250, 500, 1000, and 1500. 
For p = 0.8, values of A^ are 50, 100, 150, 200, 250, 300, 400, and 500. For p = 1.0{K = 3), 
values of A^ are 50, 100, 150, 200, and 250. 
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crossover curves for p = 0, .2, .4, .5, .6, .8, and 1 .0 
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5. Rescaled data for all p, using the largest values of N which could be obtained, with 
alphttc and u determined as described in text. 
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6. Theoretical and experimental results for the SAT/UNSAT transition in the 2+p-SAT 
model. The vertical line at po separates the continuous from the discontinuous transition. 
The full line is the replica-symmetric theory's predicted transition, believed exact for p < po, 
and the diamond data points with error bars are results of computer experiment and finite- 



size scaling. The other two lines show upper and lower bounds obtained in |2^, while the 



stronger upper bound due to JS^], and the best known lower bound, due to |3§], are indicated 
by square data points. 
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7. Crossover seen in the exponent u governing the width of the critical regime, as K 
increases from 2 to 3. 
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Fig. 8a: Median computational cost, linear scale. 



p = 0.0, 0.2, 0.4, 0.6, using tableau 
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Fig. 8b: Median cost, semilog scale. 

. Median computational cost of proving a formula SAT or UNSAT using the TABLEAU 
search method, for p ranging from to 1. The data in (a) is plotted on a linear scale, 
appropriate for the cases with p < pq. The semi-log plots in (b) show an exponential 

dependence of cost on N ior p > po. 
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Backbone fractions for k=2 ground states 
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fig. 9a: Backbone fraction for K = 2. 
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fig. 9b: Backbone fraction for K = ?>. 

9. Backbone fractions as a function of a averaged over many samples for ii' = 2 (a) cases 

with A^ = 18 to 45 and K = ?> (h) cases with A^ = 18 to 28. The vertical lines mark the 

SAT/UNSAT thresholds in the limit A^ -^ oo. For 2-SAT, data obtained from larger sizes 

A^ = 100, 200, 500 show that the backbone fraction at the threshold decreases to zero. 
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